Mechanism of Terahertz Electromagnetic Wave Emission from 
Intrinsic Josephson Junctions 



Masashi Tachiki, 1 Shouta Fukuya, 2 and Tomio Koyama 3 

1 Graduate school of Frontier Science, The University of Tokyo, 5-1-6 Kashiwanoha, Kashiwa 277-8581, Japan 
Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan 
Institute for Materials Research, Tohoku University, Katahira 2-1-1, Aoba-ku, Sendai 980-77, Japan 

(Dated: October 29, 2008) 

Using a 3-D parallelepiped model of intrinsic Josephson junctions, we calculate the cavity reso- 
nance modes of Josephson plasma waves excited by external electric currents. The electromagnetic 
(EM) wave of the excited Josephson plasma is converted to a THz EM wave at the sample surfaces. 
The cavity modes accompanied by static phase kinks of the superconducting order parameter have 
been intensively investigated. The phase kinks induce a spatial modulation of the amplitude of 
the order parameter around the kinks and decrease the superconducting condensation energy. The 
Josephson plasma produces a magnetic field in the vacuum in addition to the emitted EM wave. 
This magnetic energy detemines the orientation of the cavity mode. Taking account of the facts 
mentioned above, we obtained sharp resonance peaks in the I-V curves and sizable powers of con- 
tinuous and coherent terahertz wave emission at the cavity resonance. The emission frequencies are 
inversely proportional to the length of the shorter side of the samples in agreement with experiments. 

PACS numbers: 74.50.+r, 74.25.Gz, 85.25.Cp 



Emission of terahertz electromagnetic (EM) waves 
from intrinsic Josephson junctions (IJJ) in high tem- 
perature superconductors has been extensively studied 
EHHHHHHH- Recently, Ozyuzer et al. succeeded 
in detecting strong and continuous emission of terahertz 
EM waves from mesa-shaped samples of the high tem- 
perature superconductor Bi 2 Sr 2 CaCu208 (BSCCO) [§]. 
The general mechanism for the emission is as follows. 
When an external current is applied along the c-axis, the 
ac Josephson current in the resistive state excites a cavity 
resonance mode of Josephson plasma wave in the sample. 
The excited standing wave of Josephson plasma is con- 
verted to a terahertz EM wave at the mesa surfaces and 
the EM wave is emitted into the vacuum space. However, 
details of the mechanism have not yet been clarified, al- 
though it is important for designing the terahertz EM 
wave emitters with use of IJJ. 

Recently, X. Hu and S. Lin 1^, 11 1, and A. Koshelev 
[T3 | proposed the following new mechanism. When the 
inductive interaction between the superconducting Cu0 2 
layers in BSCCO is strong, static kink structures arise in 
the phase difference of superconducting order parameter 
between the superconducting layers. The phase kinks in- 
duce cavity resonance modes of the Josephson plasma. 
This is a new dynamic state caused by the non-linear 
effect special in the IJJ system. In this paper, we first 
discuss the stability of this new state, and then we inves- 
tigate the mechanism of the terahertz EM wave emission 
on the basis of the discussion. 

For the sample of IJJ, we use a model shown in Fig. 
[TJ In this figure the superconducting Cu0 2 layers and 
the insulating layers in the IJJ are shown in green and 
light green, respectively. An external electric current is 
applied in the direction of the z-axis, perpendicular to 
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FIG. 1: (color online) Schematic view of intrinsic Josephson 
junctions. 



the layers. The L x , L y , and L z are the sample lengths 
respectively along the x, y, and z-axes. Now, we de- 
rive the equation for the simulation. The supercon- 
ducting order parameter of the Ith layer is expressed as 
'J; = Az(r, t) exp[iipi(r, t)] with r = (x, y), x, y and t re- 
ferring to the spatial and temporal coordinates, respec- 
tively. We assume that the amplitude A;(r,t) is con- 
stant independent of space and time, and only the phase 
ipi(r,t) is dependent on space and time. In this case, 
the current density along the z-axis is given by a sum of 
the Josephson, quasiparticle, and displacement current 
densities as 

= J c sim(}i +1: i + a c E z j + n + -^-d t E z i + xi (1) 
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where J c is the critical current density, a c is the normal 
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conductivity along the c axis, and E z j+n is the electric 
field between (I + l)th and Ith layer along the z axis. In 
Eq. (1), V';+i,i( r J *) i s tne g au g e invariant phase differ- 
ence defined as 



ipi+i,i(r,t) = tpi +1 (r,t) - tpi(r,t) 



27r f Zl+1 * , 

— / A z (r,z,t)dz 



(2) 



with the vector potential A z (r,z,t) and the flux unit 
4>o ■ For the superconducting current density in the Cu02 
plane, we use the generalized London equations, since the 
Ginzburg-Landau parameter is very large in BSCCO. Wc 
insert Eqs. ([T]) and ^ into Maxwell's equations along 
with the superconducting current in the CuC>2 layers. 
Following the calculation procedure given in [6| , we have 



(l-CA (2) )(sinV/ 
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In above equation, we use normalized units for length 
and time, respectively, defined as x' = x/A c , y' = y/A c , 
t' = ujpt, where uj p is the plasma angular frequency, and 
A c is the penetration depth of the magnetic field applied 
along the x or y-axis from the xz or yz-surface. The 
parameters in Eq. ((3]) are defined as Q — \ 2 ab /sd and 
j3 = 4:ir\ c <j c /^/eZ, where the s and d are the thicknesses of 
the superconducting and insulating layers, respectively. 
X a b is the magnetic field penetration depth from the xy- 
plane surface and £ is the inductive constant between the 
CuC>2 layers. The operator A' 2 ) is defined as fi = 
fl+i—2fi+fi-i. In Eq. ©, we neglect the charging effect 
in the Cu0 2 layers, since we consider the region above the 
retrapping point, and in this region the charging effect is 
much weaker than the inductive effect between the Cu02 
layers. Keeping in mind that the IJJ is Bi 2 Sr 2 CaCu208, 
we choose A a h = 0.4^m, A c = 50 ~ 150^m, s = 3A, 
d = 121, p = 0.02, and C = 5 x 10 5 . 

We express the phase difference ipi + ij(r,t) as 



ipi+x,i(r' ,t') 



r 



(4) 



The first term in the right hand side of Eq. (HJ) is the 
phase difference due to the ac Josephson effect, the sec- 
ond term is the phase difference due to the excited cavity 
mode, the third term is a static phase difference, and the 
fourth term is the phase difference due to an external 
current density normalized by J c , I'. 

The normalized sample lengths L' x , L' y and L' z are 
taken to be 1, 3, and 0.01. Although the samples used 
in the experiments Q is composed of several hundreds of 
intrinsic Josephson junctions, to simplify the calculation 
we impose an assumption that the phase difference has 
a four junctions periodicity along the z-axis. Since the 
sample length along the z-axis is much shorter than the 
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FIG. 2: (color online) Typical configuration of the static phase 
difference. 




FIG. 3: (color online) Snapshot of the electric and magnetic 
fields at V' = 3.8 and I' = 0.09. E' z and \B'\ are electric and 
magnetic fields normalized by <f>o/2n\ c d. Fhe sample with 
the lengths L'x — 1 and L'y = 3 is used. 



wave length of terahertz EM wave, the EM waves emitted 
from the side surfaces of the sample are not plane waves. 
The boundary condition in this case has been given by N. 
Bulaevskii and A. Koshelev [lj|. We use the boundary 
condition and tentatively take B x /E z = 7 = ±0.1 at the 
xz-surfaces and B y /E z = 7 = ±0.1 at the yz-surfaces, 
E z , B x and B y being respectively the oscillating parts of 
electric and magnetic fields at the sample surfaces. Un- 
der the boundary condition, we numerically solved Eq. 
© and obtained a solution of the static phase difference 
V'f+i ii r ) with a kink and anti-kink structure at a cav- 
ity resonance voltage V = 3.8 as shown in Fig. [2] The 
normalized voltage V is defined by V/V p , V and V p be- 
ing the voltages between the CuC>2 layers and hw p /2e, 
respectively. As seen in the figure, the kink phase struc- 
ture occurs along the y-axis. The oscillating electric field 
of the cavity mode at a time is shown in Fig. [3fa) . The 
electric field E z is almost uniform along the x-axis, and 
the electric field along the y-axis is a standing wave with 
two wave lengths. The amplitude of the oscillating elec- 
tric field is symmetric with respect to the center of the 
sample along the x- and y-axes. The absolute value of the 
oscillating magnetic field |B| = JB^,+ B% at the time 

is shown in Fig. [3jb). The absolute value of B x is much 
larger than that of B y inside the sample and at the yz- 
surfaces. This state with phase kinks is similar to those 
obtained by X. Hu and S Lin 1(| U|, and A. Koshelev 
121. In the above calculation, when we took several kinds 
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of the initial condition with modulations along the z-axis 
for the phase difference ipi+ij(r, t), we obtained the sev- 
eral kinds of phase kink state. The internal energies of 
the phase kink states consist of superconducting current 
energies and electric and magnetic energies. The internal 
energies have almost the same value in independence of 
the initial conditions. If we took an initial condition that 
the phase difference is perfectly uniform along the z-axis, 
we obtained a solution of the state without phase kink 
in place of the solutions with phase kinks. The calcu- 
lated internal energies of the phase kink states are lower 
than that of the state without kink phase, and the ratio 
of the energy of the phase kink state to the energy of 
the state without phase kink is about 0.8. Therefore, the 
state with phase kinks can be excited by a lower external 
current. According to the facts mentioned above, it is 
expected that the states with phase kinks are often ex- 
cited by applying external currents. 

In the calculation so far, we have assumed that the 
amplitude of the superconducting order parameter is spa- 
tially uniform in the CuC>2 layers. However, the phase 
kinks induce a modulation of the amplitude of the order 
parameter as discussed in the following. The large change 
of the static phase difference occurs around y' = ±0.375 
and ±1.125 as seen in Fig. [2l The change causes a spatial 
change of the phase of the superconducting order param- 
eter in the superconducting CuC>2 layers. This change 
of the phase induces a reduction of the ampliude of the 
order parameter around the phase kinks. This reduc- 
tion of the order parameter amplitude is similar to that 
in the vortex core of the Abrikosov vortex. According 
to this reduction of the amplitude, the superconducting 
condensation energy is reduced and thus the energy of the 
system is increased. We estimate the energy increase in 
the following way. Using the Ginzburg-Landau equation 
we calculate the order parameter A(r) modified by the 
phase kinks. Then, the energy increase E c is calculated 
by the formula E c = {HI/Stt) f[l- (A(r) / A a ) 2 ]dV , H c 
being the thermodynamic critical field and Ao being the 
amplitude of the order parameter without phase kink. 
Using the experimental value of H c we estimate E c to be 
1^2 times as large as the internal energy of the state 
without phase kink. If E c is added to the internal energy 
of the kink phase, the total energy of the kink state is 
2^3 times larger than the internal energy of the state 
without kink. Therefore, to excite the state with phase 
kink, we need an external current larger than the current 
to excite the state without kink. If heating is considered, 
the state with kink is hardly excited [^, [lj] . 

In the calculation of the resonance mode shown in Fig. 
[3j we have missed an important fact. As shown in Fig. 
[3] the magnitude of B x is much larger than that of B y 
at the yz- surfaces. Although the magnetic field B x does 
not contributes to the radiation, it penetrates the vac- 
uum from the surface. We should add this magnetic en- 
ergy in the vacuum to the internal energy of the state. 




V 

FIG. 4: (color online) I-V curves in the two cases of (a) L' x — 2 
and (b) L' x — 3. The boundary condition parameter 7 = 0.1 
is used. 



If the resonance occurs along the x-axis instead of along 
the y-axis, the magnitude of B y becomes much larger 
than that of B x . Since the xz-surface area is narrower 
than the yz-surface area for the sample, the magnetic en- 
ergy penetrated from the rrz-surface into the vacuum is 
smaller and, thus stabilizes the magnetic field B along 
the y-axis. Therefore, the cavity resonance mode actu- 
ally occurs along the s-axis rather that along the y-axis 
in the sample [9j . The origin of the magnetic energy sta- 
bilizing B along the y-axis in this case is similar to that 
of the shape anisotopy energy of a ferromagnet that sta- 
bilizes its magnetization along the longer size direction 
of the ferromagnet. 

On the basis of these discussions, we investigate the 
physical properties originating from the cavity resonance 
along the x-axis in the states without phase kink. Figure 
0Ja) shows the calculated current- voltage (I-V) curve for 
L' x = 2. As seen in the figure, two sharp peaks appear 
respectively at the normalized voltages V = 2.9 and 6.2. 
These sharp peaks are caused by the cavity resonance 
of the excited electric and magnetic fields in the sample. 
Since the amplitude of the oscillating electric field is large 
at the resonance voltages, the sharp peaks are considered 
to be internal Shapiro steps induced by the oscillating 
electric fields in a self-consistent way [15j. Figure. HJb) 
shows the current- voltage curve for L' x = 3. The cavity 
resonance voltage at V = 1.7 almost coinsides with the 
retrapping voltage The oscillating electric and mag- 
netic fields for V = 2.9 at a time is shown in Fig. 03a). 
The amplitudes of both the electric and magnetic fields 
are asymmetric with respect to the center of the sample. 
Figure [5jb) shows the emission intensities (the time aver- 
ages of the Poynting vectors) for the boundary condition 
parameters 7 = 0.1 and 7 = 1 as functions of V'. The 
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FIG. 5: (color online) (a) Snapshot of the electric and mag- 
netic fields at the top of first internal shapiro step V' = 2.93 
and I' — 0.112. The normalization of E' z and B' y is the same 
as that in Fig. [3] (b) Voltage dependence of the intensity 
of emission from intrinsic Josephson junctions. The length of 
the sample is L' x — 2. Black circles and blue circles show the 
emission intensities in the cases of 7 = 0.1 and 7 = 1, respec- 
tively. The time averages of Poynting vectors are normalized 
by c(0 o /27rA c d) 2 /47r, 5". 
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FIG. 6: (color online) 1/L' X dependence of the first resonance 
voltage. 



black circles indicate the emission intensities for 7 = 0.1. 
The sharp emission peaks appear at the cavity resonance 
voltages as seen in Fig. E{a). The maximum emission 
power at V' = 2.9 is estimated to be 23mW for the sam- 
ple with A c = 100/im. The blue circles in Fig. [5{b) show 
the emission intensities for the boundary condition pa- 
rameter 7 = 1. The 7 = 1 means that L' y and L' z are 
much longer than the emitted EM wave length and thus 
the emitted EM waves are plane waves. As seen in the 



figure, the voltage dependence of the emission intensity is 
very different from that for 7 = 0.1 and the intensity van- 
ishes at the cavity resonance voltages. Matsumoto et al. 
have calculated the intensity in this case and shown that 
the oscillating electric and magnetic fields vanish at the 
yz-surfaces [la ]. The results indicate that the emission 
intensity strongly depends on the sample size relative to 
the EM wave length. Each of the blue circles in Fig. [6] 
denotes the lowest normalized voltage of the cavity reso- 
nance for each of 7 samples with different L' x . As seen in 
the figure, the normalized voltage V or the normalized 
cavity resonance frequency u>' is inversely proportional to 
L' x in agreement with the experiments [9j. 
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